subroutine evolve(dfu, uxMat, uzMat, Ex2D, Ez2D,  E0C, f0vec, &    ! input/output
                  vvec, pvec, zvec, kXvec, kZvec, &                ! input
                  ExFilter, phi02DFT,  fvc0T, &                    ! input
                  xMod, uPerpMod, &                                ! input
                  ExpKv, ExpKz, ExpKVdt, &                         ! input
                  E0z, Density2Dreal, &                            ! output
                  phi2DFT, &                                       ! output
                  CEz2D, EzKx,  density1DFT, &                     ! output
                  OmegaBounce )                                    ! output


  use fft_wrap, only   : fft_forward, fft_backward

  use dimensions, only : jx, jz, jv, jDirection, &   ! domain size in r-space
                         nkx, nkz, np,           &   ! limits      in k-space
                         dx,  dz,  dv,           &   ! resolution  in r-space
                         dkx, dkz, dp                ! resolution  in k-space
         

  use parameters, only : phi0, kz, vPhase,  Dxxxx, Dvvvv, &            ! real,     phys
                         ny,  iWaveFrame, iFixPump, iModx, iModz, &    ! integer,  phys
                         ModAmp,  Lx, SpeckleWidth, u, &               ! real,     phys
                         restart, nt, &                                ! integer,  num
                         Vmax, dt, KxFilter, &                         ! real,     num 
                         iself, iskip, idfMod, &                       ! integer,  num 
                         iE0Modulation, iBcPerpPeriodic, iFourier, &   ! integer,  num
                         pi, iunit, isign                              ! constants

  implicit none



  !-- input variables (from main code) --


  real,    intent(inout)   :: dfu(jz, jv, jx, jDirection)
  real,    intent(inout)   :: uxMat(jx, jz, jDirection)
  real,    intent(inout)   :: uzMat(jx, jz, jDirection)
  real,    intent(inout)   :: Ex2D(jx,jz)
  real,    intent(inout)   :: Ez2D(jx,jz)
  complex, intent(inout)   :: E0C(jz)
  real,    intent(inout)   :: f0vec(jv)
    

  real,    intent(in)      :: vvec(jv), pvec(jv), zvec(jz)
  real,    intent(in)      :: kXvec(jx), kZvec(jz)
  real,    intent(in)      :: ExFilter(jx)
  complex, intent(in)      :: phi02DFT(jx,jz)
  complex, intent(in)      :: fvc0T(jv)
  complex, intent(in)      :: xMod(jx), uPerpMod(jx)
  complex, intent(in)      :: ExpKv(jz,jv)
  complex, intent(in)      :: ExpKz(jz)
  complex, intent(in)      :: ExpKVdt
  
  real,    intent(out)     :: E0z(jx,jz)
  real,    intent(out)     :: Density2Dreal(jx,jz)
  complex, intent(out)     :: phi2DFT(jx,jz)    ! diagnostics only
  complex, intent(out)     :: CEz2D(jx,jz)      ! diagnostics only
  complex, intent(out)     :: EzKx(jx,jz)       ! diagnostics only
  complex, intent(out)     :: density1DFT(jz)

  real,    intent(out)     :: OmegaBounce 



 !-- local variables -----------------------

  complex    :: phi02D(jx,jz)

  complex    :: fz(jz), fvc(jv)
  complex    :: Phi1Dz(jz)
  complex    :: z(jv)
  complex    :: du(jx)

  real       :: aVar(jx), fVec(jx)

  real       :: fVsVvec(jv)
  real       :: Phi1D(jx)
  real       :: uTemp(jz)
  real       :: uzVec(jz), uPerpVec(jz)

  integer    :: ix, iz, iv, iDirection



 !------------------------------------------

    
   ! ExpKv(ikz,:) = cexp( -iunit*vvec(:)*kZvec(ikz)*dt - Dxxxx * kZvec(ikz)**4 *dt  )
   
   do iDirection = 1,2  
   
      do ix=1,jx
         do iv=1,jv    !	begin velocity loop
            fz(:) = dfu(:,iv,ix,iDirection)
            call fft_forward(jz,fz,fz)
            fz = fz/float(jz)
            call dealias(jz,fz) !  fz=FT(dfu)      
            fz(:) = fz(:)*ExpKv(:,iv)
            call fft_backward(jz,fz,fz)
            dfu(:,iv,ix,iDirection) = fz(:)
         enddo     !		iv
      enddo         !     ix
    
      do iz=1,jz                      
         aVar(:) = uxMat(:,iz,iDirection)           
         do iv=1,jv !  complex fvc(jv),fvc0(jv)     fvc0(iv)=f0vec(iv)
            fVec(:) = dfu(iz,iv,:,iDirection) + 0.5 * f0vec(iv)   ! real fVec(jx)
            call UpWind2ndOrder_Variable_a( fVec,  dt, dx, aVar, jx, &
                 0.5 * f0vec(iv),  0.5 * f0vec(iv), iBcPerpPeriodic )
            dfu(iz,iv,:,iDirection) = fVec(:)  - 0.5 * f0vec(iv) 
            !  0.5 f0vec  since each direction has 1/2 density,
            ! goes to 1/4 is 3D (1/6 with hexagonal symmetry)
         enddo ! iv           
      enddo       ! iz 
     
   enddo         ! iDirection
    

   call GetEfield(dfu, kXvec, kZvec, ExFilter, phi02DFT, &   ! inputs 
                  Ez2D, Ex2D, OmegaBounce, Density2Dreal, &  ! outputs
                  density1DFT, phi2DFT, &                    ! outputs
                  CEz2D, EzKx)                               ! outputs
   
     
       
     
    
   do iDirection = 1,2   !   real dfu(jz,jv,jx,jDirection)   ACCELERATION , df/dt + E df/dv = 0
      do ix=1,jx
         do iz=1,jz   !	begin spatial loop
            fvc(:) = dfu(iz,:,ix,iDirection)
            call fft_forward(jv,fvc,fvc)
            fvc = fvc/float(jv)
            call dealias(jv,fvc)
            z(:)   =  - Dvvvv*pvec(:)**4  - iunit*pvec(:)*Ez2D(ix,iz) + 1.E-09
            fvc(:) = fvc(:)*cexp(z(:)*dt) + &
                 0.5 * (1-cexp(z(:)*dt))*iunit*pvec(:)*Ez2D(ix,iz)*fvc0T(:)/z(:)
                 ! 0.5 since each direction has 1/2 density, 
                 ! goes to 1/4 is 3D (1/6 with hexagonal symmetry)
!	       do ip=1,np !  but not ip = 1 since pvec(1) = z(1) = 0 
!	         fvc(ip)=fvc(ip)*cexp(z(ip)*dt) + &
!                         (1-cexp(z(ip)*dt))*iunit*pvec(ip)*Ez2D(ix,iz)*fvc0T(ip)/z(ip)
!	       end do
            call fft_backward(jv,fvc,fvc)
            dfu(iz,:,ix,iDirection)=fvc(:)       !   f(iv, iDirection) = dfu(iv, iDirection) + 0.5 * f0vec(iv)
         end do !		end  z spatial loop
      enddo       !       end  x spatial loop
   enddo        !       end idirection loop

       
   do ix = 1,jx
      do iz = 1,jz 
         do iDirection = 1,2
            !real dfu(jz,jv,jx,jDirection), fVec(jx), fVsVvec(jv),vvec(jv), uzMat(jx,jz,jDirection)
            fVsVvec(:) = dfu(iz,:,ix,iDirection) + 0.5 * f0vec(:)
            uzMat(ix,iz,iDirection) = sum( fVsVvec * vVec ) / sum( fVsVvec )
         enddo
      enddo
   enddo
   
   
   !  ADVECT u_perp IN Z: (d/dt  - v0 d/dz) u_perp = 0 
   !  ExpKz(:) = cexp( iunit * kZvec(:) * V0 * dt )
   
   do iDirection = 1,2
   
   
      IF(iWaveFrame==1)  THEN  !  advect ux in z since <u_z> approximately = -vphase =  V0, added on June 2, 2009
         do ix=1,jx  !  complex Phi1Dz(jz), ExpKz(:) = cexp( iunit * kZvec(:) * V0 * dt )
            if( iFourier==1 ) then    !  Fourier with linearization in uz

               Phi1Dz(:) = uxMat(ix,:,iDirection)  !  real uxMat(jx,jz,jDirection)
               call fft_forward(jz,Phi1Dz,Phi1Dz)
               Phi1Dz=Phi1Dz/float(jz)
               call dealias(jz,Phi1Dz) 
               Phi1Dz(:) = ExpKz(:) * Phi1Dz(:)
               call fft_backward(jz,Phi1Dz,Phi1Dz)
               uxMat(ix,:,iDirection) = Phi1Dz(:)
            endif
            if( iFourier.EQ.-1 ) then   !  upwind
               uzVec(:) =    uzMat(ix,:,iDirection)  !  advecting velocity in z direction, at each x
               uPerpVec(:) = uxMat(ix,:,iDirection)
               call SubUpwind_NonCon_Form( uPerpVec, dt, dz, uzVec, jz, 0.0,  0.0, 1, 1, uTemp ) !  periodic in z so far
               uxMat(ix,:,iDirection)= uPerpVec(:)
            endif
            if( iFourier.EQ.0 ) then   !  centered difference
               uzVec(:) =    uzMat(ix,:,iDirection)  !  advecting velocity in z direction, at each x
               uPerpVec(:) = uxMat(ix,:,iDirection)
               call SubUpwind_NonCon_Form( uPerpVec, dt, dz, uzVec, jz, 0.0,  0.0, 1, 0, uTemp ) !  periodic in z so far
               uxMat(ix,:,iDirection)= uPerpVec(:)
            endif
         enddo   ! real  uTemp(jz), uzVec(jz), uzMat(jx,jz,jDirection), , uPerpVec(jz)
      ENDIF  ! else need to impliment centered difference variant of SubUpwind_NonCon_Form
     
      uxMat(:,:,iDirection) = uxMat(:,:,iDirection) + Ex2D(:,:) * dt  ! uxMat(ix,iz,iDirection)  all ix and all iz
         
   enddo        !       end idirection loop
 


   do iDirection = 1,2     !  transverse self-advection of u_perp       data iSign(1), iSign(2)/ 1, -1 /
   
      do iz=1,jz  !  real f0vec(jv)
         fVec(:) = uxMat(:,iz,iDirection)  !  real fVec(jx), uxMat(jx,jz,jDirection)
         call NonLinearAdvection( fVec, dt, dx, jx, u*iSign(iDirection),  u*iSign(iDirection), iBcPerpPeriodic )  
         !  MAYBE NOT periodic bc for ux
                      
         du = fVec - u * iSign(iDirection) !  approx to uniform part of du,  improve upon this         complex du(jx)
         du = du * uPerpMod  !  make du periodic in x
         call fft_forward(jx,du,du)
         du=du/float(jx)
         call deAlias(jx,du)
         du = du * ExFilter
         call fft_backward(jx,du,du)
         fVec = du +  u * iSign(iDirection)                      
         uxMat(:,iz,iDirection) = fVec(:)
      enddo  ! iz
          
   enddo      ! iDirection loop                     data iSign(1), iSign(2)/ 1, -1 /
                              
     
         
   E0C(:) = ExpKVdt * ExpKVdt * E0C(:)
         
   if (iE0Modulation==1) then               !  modulate  source in x direction
      do ix=1,jx   
         ! E0C(iz) =  -iunit * kz * phi0 *exp ( iunit * kz * (iz-1) * dz  )
         E0z(ix,:) =      xMod(ix) * E0C(:) 
         ! so when x=SpeckleWidth/2, E0z is down by a factor of 0.00674
         phi02D(ix,:) =   xMod(ix) * phi0 *exp ( iunit * kz * zvec(:)  ) 
         ! E0C(iz) = -iunit * kz * phi0 *exp ( iunit * kz * z  )
      enddo
   else
      do ix=1,jx
         E0z(ix,:) = E0C(:)
         phi02D(ix,:) =  phi0 *exp ( iunit * kz * zvec(:)  )
      enddo
   endif
    
    

end subroutine
